pub fn find_max_prime_factor(number: u128) -> u128 {
    if number <= 1 {
        return number;
    }

    let mut n = number;
    let mut max_factor = 0;

    while n % 2 == 0 {
        max_factor = 2;
        n /= 2;
    }

    // For small primes, quickly check and divide
    let small_primes = [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47];
    for &p in &small_primes {
        while n % p == 0 {
            max_factor = p;
            n /= p;
        }
    }

    // Fast modular exponentiation
    fn mod_pow(mut base: u128, mut exp: u128, modulus: u128) -> u128 {
        if modulus == 1 {
            return 0;
        }
        let mut result = 1;
        base %= modulus;
        while exp > 0 {
            if exp & 1 == 1 {
                result = result * base % modulus;
            }
            base = base * base % modulus;
            exp >>= 1;
        }
        result
    }

    // Miller-Rabin primality test
    fn is_prime(n: u128) -> bool {
        if n <= 1 || n == 4 {
            return false;
        }
        if n <= 3 {
            return true;
        }
        if n % 2 == 0 {
            return false;
        }

        let bases = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37];
        let mut d = n - 1;
        let mut r = 0;
        while d % 2 == 0 {
            d /= 2;
            r += 1;
        }

        'next_base: for &a in bases.iter().take_while(|&&x| x < n) {
            let mut x = mod_pow(a, d, n);
            if x == 1 || x == n - 1 {
                continue;
            }
            for _ in 0..r - 1 {
                x = x * x % n;
                if x == n - 1 {
                    continue 'next_base;
                }
            }
            return false;
        }
        true
    }

    // Quadratic probing
    fn find_factor(n: u128) -> u128 {
        if n % 2 == 0 {
            return 2;
        }
        if is_prime(n) {
            return n;
        }

        let mut x = 2u128;
        let mut y = 2u128;
        let mut d = 1u128;
        let mut i = 0;

        while d == 1 {
            i += 1;
            x = (x * x + i) % n;
            y = (y * y + i) % n;
            y = (y * y + i) % n;

            if x > y {
                d = gcd(x - y, n);
            } else {
                d = gcd(y - x, n);
            }

            if d == n {
                x = (x + 1) % n;
                y = (y + 1) % n;
                d = 1;
            }
        }
        d
    }

    fn gcd(mut a: u128, mut b: u128) -> u128 {
        while b != 0 {
            let t = b;
            b = a % b;
            a = t;
        }
        a
    }

    while n > 1 {
        if is_prime(n) {
            max_factor = max_factor.max(n);
            break;
        }
        let factor = find_factor(n);
        max_factor = max_factor.max(factor);
        n /= factor;
    }

    max_factor
}
